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Abstract 

We develop two further numerical tests of parton cascade models: the 
ideal gas equation of state and the collision frequency tests. The equa- 
tion of state test checks the initial momentum distribution generator 
and free expansion dynamics in periodic inhomogeneous geometries. 
The collision rate test is sensitive to spatial inhomogeneities and the 
collision algorithm. These tests are applied to the recently developed 
ZPC parton cascade model. The tests helped uncover unphysical cor- 
relations induced by one of the commonly used random number gener- 
ators and showed the necessity of particle subdivisions for convergence 
to the exact analytic limit. 
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1 Introduction 



Many observables in nuclear collisions are difficult to calculate analytically 
because the number of particles is neither large enough to justify rigorously 
the application of statistical mechanics nor small enough to justify impulse 
approximations. In most cases, numerical simulations of transport equations 
are required to compare theory with experiment. However, due to the com- 
plexity of the algorithms employed and the number of untested dynamical 
assumptions 0, |3|, Q , it is not straight forward to check or even reproduce 
the numerical results from parton cascade codes. Recently, several analytic 
tests were proposed to help check the accuracy and identify limitations of 
such codes. In the context of earlier non-relativistic transport models tests 
such as comparing with the analytic Krook-Wu model || have been use- 
ful. In the context of hydrodynamics, numerical hydrodynamic codes have 
been tested in cases of expansion of baryon free matter into vacuum || and 
"slab-on-slab" collision j7]. 

To address the new generation of parton cascade, a new Open Standard 
for Codes and Routines (OSCAR) || has been developed to enable objective 
testing of essential components of algorithms and ensure reproducibility of 
numerical results. In this context, we have proposed in ref.0| a dynamical 
test of the evolution of the transverse energy in 1 + 1 dimensional expansion 
for cascade models. Another test of the frame and scattering scheme depen- 



dence of cascade models was proposed in [|T(], [nfl . In this paper, we propose 
two further tests for the OSCAR standard: the equation of state and the 
collision rate test of one-component relativistic classical gases. 

The tests provide information about the nature of the code's initial mo- 
mentum distribution and check its evolution algorithm in free expansion in 
periodic spatial inhomogeneities against nontrivial analytic expressions. The 
transverse to parallel momentum flux ratio tests the evolution of space-time 
and momentum space correlations. The asymptotic homogeneous equili- 
brated state tests the equation of state of the model against the expected 
ideal gas laws. The collision frequency test checks the basic collision algo- 
rithm. 

We apply these tests to the newly developed ZPC parton cascade model 
[|]]. In order to pass these tests, we uncovered and fixed a problem with one 
of the random number generators used in ZPC. We discuss that example 
in detail to illustrate the importance of applying such numerical tests. We 
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propose to add these tests to the group already in the OSCAR standard |8|]. 

The paper is organized as follows: In section 2, we recall basic statistical 
mechanics relations of ideal gas and calculate the free streaming evolution 
of periodic slab initial spatial distribution. In addition the collision rate is 
discussed. Numerical tests of ZPC are compared to analytical predictions in 
section 3. We conclude by emphasizing the importance of testing cascade 
codes. 



2 Analytic Tests 
2.1 Equation of State 

The partition function for an ideal relativistic classical gas with mass m and 
degeneracy 7 is given by J15] (h = c = k% — 1) 

Z(N,T 1 V) = ±(^(Tm 2 )K 2 (m/T)e m / T ) N , (1) 



AH V2tt 2 

In the N ^> 1 limit, the free energy, F = —TlogZ is well approximated by 

F(N,T,V) = -N\T\og(^^ r K 2 (m/T))+m + T\ . (2) 

The pressure, P = —dF/dV is given by the well known ideal gas law 

P(N,T,V)=pT , (3) 

in terms of the density p = N/V, and the energy density is given by 

e(N, T, V) = pm (K^m) / K 2 ((3m) + 3T/m) . (4) 

where (3 = 1/T. For fixed N, V an interesting quantity that reflects the 
softness of the equation of state is 

- = 1 (5) 

e 3 + PmK x {^m)/K 2 {(3m) ' v ' 

In the m/T <ti 1 relativistic limit, P/e w 1/3, while in the nonrelativistic 
limit, P/e w T/m, 
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In a cascade simulation, these thermodynamic quantities can be measured 
as a function of time by computing the spatially averaged energy-momentum 
tensor 

<T^>=W^- (6) 

Then e = (T 00 ) and P = (T 11 ). In section 3, we will compare the analytic 
results of the pressure-energy density ratio as a function of m/T, and the 
energy density as a function of T. 



2.2 Free Expansion In Slab Geometries 

A simple test of equilibration in a periodic box of volume V = L 3 and fixed 
particle number is provided by taking the initial spatial distribution to be 
confined to one half of the box (say with x > 0): 

f(x,p,0) = Po (l + -Imtanh-^e 2 ^])/^). (7) 

71 

In the above, po = N/V is the equilibrium particle density, L is the length 
of the box, and 

, exp(-/Vm 2 + p 2 ) 
J[P) 47r(m 2 / '(5)K 2 {(3mY {) 

For the free streaming case, f(x,p,t) = f(x — pt/p°,p,0). 

For the interesting case of massless partons, the momentum flux compo- 
nent, T 11 evolves from its initial value 

T u (x, 0) = T n (oo)(l + -Imtanh-^e 2 ^))] (9) 

7T 

to its final value T n (oo) — po I d 3 p pl/p°f(p) via 

^4 = 1 + Ts £ ^%^(((^) 2 - 2 ) sin (^) + 2 tocos(M)- (10) 
T (°°) * ^ n 

Here a = 3L 3 /(27r 4 ), n = 2 / + 1, and fc n = 2rm/L. The infinite sum can be 
expressed in terms of the Lerch n, a) functions, but in practice the series 
converges very rapidly. For example the evolution at the "midway" points 
x = (2k + l)L/4 involves damped oscillations with maxima and minima 
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occuring at times t = (2m + l)L/4. At x — L/4, the "11" momentum flux 
oscillates within the envelops 



r»(, = L/4,0)=T»( OO )(l±i("-^)) I") 

reaching the envelop at times t = (2m + l)L/4. Note that free streaming in 
this geometry leads to global homogeneous equilibrium through mixing from 
neighboring slabs. The deviation of the parton cascade evolution from this 
result when collisions terms are turned on are of physical interest as tests of 
collective hydrodynamic behavior. 

Another interesting probe of the evolution is the transverse-longitudinal 
momentum flux anisotropy, 

A zx (x,t) = (T 33 (x,t))/{T u (x,t)) (12) 

Initially in the — L/2 < x < region we set A zx = while in the < x < L/2 
region A zx = 1. In the case of free streaming, this anisotropy is given by 

AZX(x f) 1 + f Eg ^(sin(U)-Ucos(A; TO t)) (u) 

1 + f E£ £ ^(((M 2 -2)sin(U)+2A; n tcos(A; n t))' 

We note that of course other non-equilibrium initial conditions, e.g., ho- 
mogeneous spatial distribution with anisotropic momentum distribution, can 
be set up to test dynamical relaxation toward global equilibrium. In this pa- 
per we will study only the evolution of spatial inhomogeneous conditions. 

2.3 Collision rate for an isotropic, homogeneous gas 

The collision rate per unit volume W is another quantity that can be easily 
monitored. It is related to the number of collisions N c in a period of time t 
via: 

w = w (14) 

We take the integrated parton elastic cross section to be: 

°99—99 = ~2- (15) 
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The total collision rate per unit volume can be calculated from: 



W = < apip 2 \v\ - v 2 \ > 

1 7T , f d 3 Pl d 3 P2 0(El+E2)] 



2^ 7 J (2vr) 3 (27r) 3 11 2| ' K > 

and it is expected to be independent of the microscopic form of the differen- 
tial cross section. The 1/2 on the right hand side of equation (|T6|) takes into 
account the two identical incoming particles. The identity of two incoming 
gluons has been taken into account in the scattering cross section. In ap- 
pendix A, we show that the above 6-dimensional integral can be reduced to 
a 1-dimensional integral: 

TT jj, 1 

in which F(x) = dy y 2 (y 2 - x 2 )K x {y). 



3 Numerical Tests of the ZPC model 

For the cascade calculations, we take the temperature T = 0.5, 1, 1.5 GeV, 
and set the total number of particles in the periodic box to be iV = 4000. 
The volume is related to the density of particle at zero-chemical potential, 

"=(2^ 4 'f dppV " E (18) 

via V = N/p. In ZPC, there are 1000 computation cells of volume V/1000. 
We specify the cell size to the cascade code. 

To study the mass dependence, we choose 3 different masses: 0, 2.5, and 
5 GeV. For the dependence on the scattering cross section, the interaction 



length ya/n = 1/pis set to be 0.5A, A, and 1.5A. In the above A is a rough 
estimate for the mean free path, 

A = — . (19) 
pa 

We specify the screening mass p to the cascade code. 
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The parameters for the cascade measurement are shown in Table (|T]). 

Fig. 1 shows the pressure-energy density ratio P/e as a function of m/T, 
and Fig. 2 shows the energy density as a function of temperature T for dif- 
ferent particle masses. We see very good agreement between the predictions 
and the cascade results. This indicates that the initial momentum distribu- 
tions are correctly generated. There is no time dependence of the pressure 
and energy density over a time period of 6 fm during which each particle has 
experienced ~ 10 collisions on average. 

For periodic slab initial conditions, T 11 evolution at position L/4 in the 
free streaming case is compared with prediction in Fig. 3. Also shown in 
Fig. 3 is the result (filled circles) with interactions turned on. We see that 
interactions reduce the rate of collective expansion compared to free stream- 
ing and also the oscillations damp faster than in the free streaming case. 
For the non-interacting free streaming case, good agreement of prediction 
(dashed line) and cascade results (pluses) is shown for the T 11 evolution at 
L/4 in Fig. 3 and time evolution of free streaming T 11 spatial distribution 
in Fig. 4. Time evolution of free streaming A zx spatial distribution (Fig. 5) 
also agrees well with the prediction. 

Fig. 6, 7 give scaled collision rate per unit volume as a function of m/T. 
The data points with the same m/T, and same \x overlap. The three data 
sets at the same m/T correspond to three different screening masses \i. The 
smaller the screening mass, the lower the collision rate. 

The collision rate is the same for the ZPC default Yukawa type scattering 
differential cross section and the straight line propagation. This indicates the 
collision number loss is not due to particle shielding. 

With small screening mass, the interaction range is large. When the in- 
teraction range is much larger than the mean free path, non-causal collisions 
become more abundant. To process the non-causal collision, we pick up one 
collision out of several collisions for one particle according to the ordering 
time (see ref. (TTJ for details). This process neglects other collisions in the 



same non-causal collision set. Some of them will not be recovered later. The 
larger the interaction range, the larger percentage of non-causal collisions. 
Hence, more collisions are neglected and the collision rate is lower than ex- 
pected. 

In the dilute limit, the percentage of non-causal collisions out of total 
number of collisions for massless particles is proportional to the number of 
particles inside the causal sphere in the two colliding particle center of mass 
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frame. The radius of the sphere is proportional to the impact parameter. So 



N nonc 1 /VI „ UA Air u3 _ 8 _ a 3 / 2 
oc — / znbdb — o 7p = — 7/?- 



iVtotal cr Jo 3 15 V7T 

2ixbdb/a is the probability of having impact parameter b. Here, 7 is for the 
averaged boost factor from lab frame to the two colliding particle center of 
mass frame. In the case of massless particles, 7 ~ 2. 47r6 3 7p/3 gives the 
number of particles in the sphere with radius b in the two colliding particle 
center of mass frame. The exact radius of the sphere depends on the definition 
of non-causal collisions [l], [TT[] and the collision prescription of cascade. 



The non-causal to total ratio is closely related to the ratio of interaction 
length to the mean free path: 

^3/2 



X = -T- = P 



pa 



We see iV nonc /iV tota i decreases linearly with x as x —> 0. This motivates the 
algorithm of reducing the non-causal collision percentage [JIJ by subdividing 
the particles by a factor I so that p — > Ip while decreasing — > a /I. This 
preserves the mean free path 1/crp while xccl/v^^O. 

The number of non-causal collisions, and total number of collisions for 
4000 particles with m = 0, T = 0.5 GeV, and a time period of 6 fm is 
summarized in table (|2]). In the table, /1 is screening mass, a gives the ratio 
of interaction length to the estimate of mean free path before rescaling. I 
is the scaling parameter, i.e., the total number of particles is increased by a 
factor of I and the cross section is decreased by a factor of I (but the number 
of collisions is still for 4000 particles). Xi is the ratio of interaction length 
to the estimate of mean free path including the rescaling, i.e., Xi — 
X2 is the ratio of interaction length to the measured mean free path. The 
mean free path is measured through the formula, p/(2W), in which W is the 
collision rate per unit volume and the particles are moving at the speed of 
light. The ratio of the number of non-causal collisions to the total number 
of collisions is also plotted in Fig. 8 against the ratio of interaction length to 
the mean free path. In Fig. 8, the open circles are data against Xi an d the 
filled circles are against X2- It shows that when the density increases, i.e., 
when x increases, the difference between X2 and Xi increases. This tells us 
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that when density is high, the naive formula for the estimate of mean free 
path, A = l/(pa), needs to be corrected. Also we see the data deviate from 
linear formula and has a tendency of saturation when density is large. This 
is consistent with the fact that the ratio should always be less than 50% from 
the definition of the non-causal collision. 

When we fix n and increase / from 1 to 5, the total number of collisions 
goes up from 64700 for / = 1 to 75400 for / = 5. For I — 10 , the total number 
of collisions is 77800. We see clearly the trend toward a constant value of 
total number of collisions when I is increased. The collision rate with I = 10 
is shown in Fig. 6 as an open triangle (the / = 1 data is shown as an open 
circle). The collision rate with / = 10 is within 1% of the analytic prediction. 

During the preliminary study of the collision rate, we found when m/T = 
10, the collision rate is higher than the predicted rate. By looking more 
carefully into the code, we found out that it was caused by the larger than 
statistical fluctuations in the position distribution. When the fugacity A is 
not uniform in space, the collision rate per unit volume: 

W = j ^\\x)I, (20) 

in which: 

By using the inequality, 

jP^>_{lp (x) ) 2 , (22 , 

and the fact that the system we prepared has zero chemical potential and 
hence averaged fugacity is one, we arrive at: 

The equality holds only when the fugacity has no spatial dependence. This 
shows that when there are space clusters existing for some time period, the 
collision rate is higher than that expected for the uniform system. 
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It was found only in m/T = 10 and not in other cases because in the 
m/T = 10 case, particles are moving very slowly, and they stay in clusters 
for much longer time. 

We traced the origin of the non-uniform distribution. It was caused by 
some correlation of random number generators (see Appendix B). When we 
use rani from |17 |, and generate first the momenta for all the particles and 
then generate the positions for all the particles, there are no abnormal fluc- 
tuations. When we generate momentum and position together, we found 
abnormal fluctuations. This does not occur when ran3 from jLTj] is used. We 
correct the generation by separating the generation of particle momentum 
and particle position. 



4 Conclusions 

From the above study, we show that the equation of state and the collision 
rate can be used to test the initial conditions and collision mechanisms of 
relativistic parton cascade. For massless particles, when the interaction range 
is much larger than the mean free path, the cascade collision rate is lower 
than the theoretical value. Other methods, e.g., particle partition, have to 
be used to correct the collision rate. 

The comparison of free streaming and interacting cascade approach to 
equilibrium indicate qualitative similarities of the two cases. However, the 
damping and speed of collective motion are quite different. A detailed com- 
parison of free streaming, ideal hydrodynamics, and cascade approach to 
global equilibrium in the case of half filled periodic box initial conditions will 
be addresses in another paper [0 . 



As discussed in this paper, spatial distribution with larger than statistical 



fluctuations gives higher than thermal reaction rate. HIJING [13| predicts 



initial spatial clusters of partons for nucleus-nucleus collisions at collider 



energies. This implies higher than thermal collision rates |14]] and many other 
interesting physical phenomena beyond the widely used hot gluon scenario 
[ |16|1 predictions. 

We emphasize the importance of using analytic tests in debugging numer- 
ical simulation codes. The spatial distribution with abnormal fluctuations 
illustrates well the usefulness of the analytic collision rate test. As more 
components are added to the cascade code, more tests will be needed to 
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ensure the consistency of different parts of the cascade code and to enable 
disentangling of the actual physical assumptions that define the model. 
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A Reduction of the phase space integral in 
the collision rate calculation 

To calculate the collision rate per unit volume, equation (|T6|), we first replace 
\vi -v 2 \ by 



s(s — Am 2 ) 
2E\E 2 

and get: 

W = --7 2 r ds Js(s-Am 2 ) f-*?* 

2 /i 27 i 4m2 V 1 ' J (2tt) 8 2£i (2tt) 3 2£ 2 

e-^ +E ^S(s-( Pl + P2 ) 2 ). (24) 
The delta function can be used to integrate out one angle by using: 

(25) 



5(s- (pi+p 2 y 



I f s-2m 2 -2E 1 E 2 
-0 I cosy 2 H — 



2pip 2 V " ^PlP2 ) 

After carrying out the angular integrals, the collision rate becomes 

-Airpl dpi 2i^p\ dp 2 



W = 



7T7 



4(2tt) 6 /i 2 J 4m' 



/ ■ 

ds ys(s — Am 2 )- 



E X E 2 2p x p 2 



i+^) ( ! _ 



s - 2m 2 - 2E l E 2 
2piP2 



(26) 
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The G function constraint can be written as: 

s 2 

{E l + E 2 ) 2 m 2 < sE x E 2 + m 2 s - —. 

Now we change integration variables from E±, E 2 to x = E\E 2 , and 
y = Ei + E 2 . The collision rate becomes: 



W 



7cy 2 f°° j— — — f e y l T 2 s y 2 m 2 , 



r°° i r e — yi 

/ ds \/s(s-4m 2 ) / dxd-u . = =Q(x + m 2 
hm? v J \/v 2 - Ax 



2(2tt) 4 /i 2 -W V v 'J » ^jyi _ Ax v 4 s ' 

(27) 

The integration over x can be carried out first. The absolute lower bound 
of x = E\E 2 is m 2 , but x is also restricted by the function. By noticing 
y > s, and hence y 2 m 2 /s — m 2 + s/4 > m 2 , we see the lower bound should 
be y 2 m 2 / s — m 2 + s/4. The upper bound is determined from the square root 
in the integrand to be y 2 /A. The result is: 

/•S/ 2 /4 1 1 i 

/ dx . . = = - J(s-Am 2 ){y 2 /s-l). (28) 

Now the collision rate per unit volume is reduced to a 2-dimensional 
integral: 



W ^ 



roc roo > 

/ ds (s - 4m 2 ) / dy e- y/T Jy 2 - s. (29) 



4(277) V 2 J4m2 JV~s 

The y integral can be readily carried out to be: 



/J dy e-y' T ^l = y/sTKx ^ . 



(30) 



In the above, K\(x) is the modified Bessel function. 
Taking 7 = 16, we arrive at equation dl7f) . 



B One example showing the interference of 
generations of distributions of variables 

The following is the a sample program that generates the particle momentum 



and position together, rani is taken from [O and declared as real*8 function 
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program dist 

implicit real*8 (a-h, o-z) 

parameter (mul = 40000) 

parameter (size = 5 . 19653401d0) 

external rani 

common /para/ xmp 

temp = 0.5d0 

xmp = 5d0 

do i = 1 , mul 

call energy (e, temp) 

call momentum (px, py, pz, e) 



X 


= 2d0 


* ranl(iseed) 


- IdO 


X 


= x * 


5d0 * size 




y 


= 2d0 


* ranl(iseed) 


- IdO 


y 


= y * 


5d0 * size 




z 


= 2d0 


* ranl(iseed) 


- IdO 


z 


= z * 


5d0 * size 





end do 

stop 

end 



subroutine energy (e, temp) 
implicit real*8 (a-h, o-z) 
external rani 
common /para/ xmp 
1000 continue 

e = ranl(iseed) 
e = e * ranl(iseed) 
e = e * ranl(iseed) 
if (e .le. OdO) goto 1000 
e = - temp * log(e) 
if (ranl(iseed) .gt. 
& exp((e - sqrt(e ** 2 + xmp ** 2)) /temp)) then 
goto 1000 
end if 
return 
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end 



subroutine momentum (px, py, pz, e) 
implicit real*8 (a-h,o-z) 
external rani 

parameter (pi = 3 . 14159265358979d0) 

cost = 2d0 * ranl(iseed) - IdO 

sint = sqrt(ldO - cost ** 2) 

phi = 2d0 * pi * ranl(iseed) 

px = e * sint * cos (phi) 

py = e * sint * sin(phi) 

pz = e * cost 

return 

end 



The momentum distribution is shown in Fig. 9 along with that generated 
by ran3, and the theory prediction. We see that the momentum distribution 
is correctly generated with reasonable fluctuations. 

Fig. 10 gives the position x distribution from the above program using 
rani and ran3. rani result has much larger fluctuations than the expected 
statistical fluctuations, while ran3 result is consistent with the expected sta- 
tistical fluctuations. 

When the generation of momentum and position are separated, the dis- 
tributions all have reasonable fluctuations. Another possible (and more ef- 
ficient) way to solve the problem is to generate Gaussian momentum distri- 
bution for large m/T values. 

Since with the above program, the total number of random number used 
is 141931896, while we can easily show the period of rani with the given 
parameters is beyond 2000000000, it is not clear that the large fluctuations 
are due to the period of the random number generator rani. But this is a 
concrete example of ill-generated distribution. 
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Table 1: Input parameters to ZPC: cell size and screening mass for differ- 
ent temperatures and masses of particles. Hi: H2, A*3 correspond to three 
decreasing screening masses as described in the text. 



T-ni; OeV 

-L W.tJ V.I C V 


m — D GeV 

111 W VI I V 


m — 2 5 GeV 

111 Zj . xj VJfC V 


m — 5 GeV 

111 tj . V.I c~ v 


Ic (fm) 
A*i (fm -1 ) 
Hi (fm -1 ) 
/i 3 (fm -1 ) 


0.53240597 
5.50178669 
4.36677096 
3.46590838 


1.31503518 
2.22745682 
1.76793364 
1.40320985 


5.19653401 
0.56368034 
0.44739338 
0.35509636 


T = 1.0 GeV 


m = GeV 


m = 2.5 GeV 


m = 5. GeV 


k (fm) 
Hi (fm -1 ) 
/i 2 (fm -1 ) 
^3 (fin -1 ) 


0.26620298 
11.00357361 
8.73354210 
6.93181691 


0.36766446 
7.96700354 
6.32341486 
5.01889767 


0.65751757 
4.45491373 
3.53586735 
2.80641975 


T = 1.5 GeV 


m = GeV 


m = 2.5 GeV 


m = 5. GeV 


Ic (fm) 
/ii (fm -1 ) 
/i 2 (fm -1 ) 
H 3 (fm -1 ) 


0.17746865 
16.50536062 
13.10031331 
10.39772550 


0.20981078 
13.96107515 
11.08091261 
8.79492611 


0.29320839 
9.99011002 
7.92915553 
6.29337487 



Table 2: Measurement of non-causal collision number and total collision 
collision number. 



H (fm : ) 


a 


1 


Xi 


X2 


-Vnonc 


A'total 


R 


5.50178669 


0.5 


1 


0.5 


0.48 


3240 


31400 


10.3% 


4.36677096 


1 


1 


1 


0.88 


8250 


46100 


17.9% 


3.46590838 


2 


1 


2 


1.56 


19000 


64700 


29.4% 


3.46590838 


2 


5 


0.89 


0.81 


12500 


75400 


16.6% 


3.46590838 


2 


10 


0.63 


0.59 


9700 


77800 


12.5% 
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Figure 1: Pressure-energy density ratio P/e as function of m/T. The dotted 
curve is the analytic prediction, and the filled black dots are cascade data. 
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Figure 2: Energy density e as a function of temperature T for different 
The curves are predictions, and the dots are cascade data. 
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Figure 3: Evolution of the longitudinal momentum flux. T 11 is in units of 
T n (oo). The size of the box is determined by the asymptotic temperature 
T = 1.5 GeV, and the screening mass is set to be 16.5 fm -1 for the interacting 
case. 
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Figure 4: Evolution of T 11 spatial distribution. The predictions are drawn 
with solid and dashed lines alternatively with numbers indicating the times 
in fm. The cascade data are drawn in filled circles and pluses alternatively 
for different times. Asymptotic temperature T = 1.5 GeV. 
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Figure 5: Evolution of A zx spatial distribution. The predictions are drawn 
with solid and dashed lines alternatively with numbers indicating the times 
in fm. The cascade data are drawn in filled circles and pluses alternatively 
for different times. Asymptotic temperature T = 1.5 GeV. 
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Figure 6: Collision rate per unit volume W as a function of m/T . The curve 
is prediction, and the pluses are cascade data. The open circle corresponds 
to the rate at m/T = with /x = 3.46590838 fm -1 . The result from scaling 
the cross section down by a factor of 10 and increasing the particle density 
by a factor of 10 is shown as the open triangle. 
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Figure 8: The number of non-causal collisions over the total number of col- 
lisions as a function of the ratio of interaction length to the mean free path. 
The dashed curve is the estimate when the radius of causal sphere is taken to 
be the interaction range. Open(filled) circles are data plotted against Xi(X2) 
in table @). 
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Figure 9: Momentum distribution for m = 5 GeV, T = 0.5 GeV with 40000 
particles. The dotted line is the prediction, the open circles come from code 
using rani, and the filled circles are ran3 results. 
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Figure 10: Number of particles per bin as a function of position x for m = 5 
GeV, T = 0.5 GeV with 40000 particles . The dotted curve is rani result, 
while the solid curve is for ran3. 
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